Quantum annealing of the Traveling Salesman Problem 
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We propose a path-integral Monte Carlo quantum annealing scheme for the symmetric Traveling 
Salesman Problem, based on a highly constrained Ising-like representation, and we compare its 
performance against standard thermal Simulated Annealing. The Monte Carlo moves implemented 
are standard, and consist in restructuring a tour by exchanging two links (2-opt moves). The 
quantum annealing scheme, even with a drastically simple form of kinetic energy, appears definitely 
superior to the classical one, when tested on a 1002 city instance of the standard TSPLIB. 
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Quantum annealing (QA) - that is using quantum me- 
chanics to optimize, through annealing, hard problems 
of everyday life, including those that have nothing to do 
with quantum mechanics - is a relatively new and fas- 
cinating idea, with important physical implications and 
potential impact in a variety of areas, from technologi- 
cal applications to other disciplines of science, wherever 
optimization of a complex system is the issue. 

The idea of QA is an offspring of the celebrated ther- 
mal simulated annealing (SA) ^ |3| , where the problem 
of minimizing a certain cost (or energy) function in a 
large space of configurations is tackled, through a sta- 
tistical mechanics analogy, by the introduction of an ar- 
tificial temperature variable which is slowly reduced to 
zero in the course of a Monte Carlo (MC) or Molecular 
Dynamics simulation. This device allows to explore the 
configuration space avoiding trapping into local minima, 
often providing a more effective and less biased search 
for the minimal "energy" than standard gradient-based 
minimization methods. 

Why not using quantum mechanics for the same pur- 
pose? Quantum mechanics works with wave-functions 
that can equally well sample wide regions of phase-space. 
Instead of thermal fluctuations, one exploits here the 
quantum fluctuations provided by a suitably introduced 
- and equally artificial - kinetic energy. Annealing is 
then performed by slowly reducing to zero the amount 
of quantum fluctuations introduced. Quantum fluctua- 
tions have, in many respects, an effect similar to that of 
thermal fluctuations - they cause, for instance, solid he- 
lium to melt even at the lowest temperatures - but they 
differ considerably in other respects. In particular, quan- 
tum systems can tunnel through classically impenetrable 
potential barriers between energy valleys, a process that 
might prove more effective than waiting for them to be 
overcome thermally, as in SA. 

Formulated in the mid 90's 0, the idea of QA picked 
up momentum only recently, through experimental work 
such as that of Brooke et al. ^ |^ , where it was shown 



that in the disordered Ising ferromagnet LiIIoo.44Yo.56F4, 
QA is both experimentally feasible and apparently su- 
perior to thermal annealing. Stimulated by these find- 
ings we recently carried out a benchmark comparison of 
classical and path-integral MC annealings on the two- 
dimensional random Ising model 0, 13- In that study, 
we confirmed the superiority of QA. We also presented 
a simple theory based on Landau-Zener tunneling for 
that. While other theoretical efforts had also reported 
some success with Q A for other problems |^ 0, ^3 ^3 > 
it is nonetheless fair to say that our overall experience 
with tackling hard problems by QA is still very limited. 
The Traveling Salesman Problem (TSP), a classic hard 
optimization problem, provides an ideal playground for a 
further test of QA in comparison with SA. In this Letter 
we report an application of QA to TSP, where we find 
that again it is superior to SA. 

Given N cities with set inter-city distances dij, TSP 
consists in finding the shortest route connecting them, 
visiting each city once and returning to the starting point. 
The literature on TSP is vast, and e.g. Ref. can be 
consulted for an account of the algorithms proposed. SA, 
although never a winning scheme when compared to some 
of the ad-hoc algorithms specifically tailored to the TSP 
problem, is known to be a very fiexible, simple, and com- 
petitive algorithm for the TSP as well The natural 
question is now: can QA do better than SA, with a com- 
parable computational effort? We will now show that 
this is indeed the case. 

The crucial point in devising a QA scheme is how to 
describe the Hilbert space of the problem, and how to 
write and implement a quantum Hamiltonian Htsp — 
Hpot + Hkin ■ Here Hpot represents the classical potential 
energy of a given configuration (in our case, the length 
of a tour), and Hkm is a suitable kinetic energy operator 
providing the necessary quantum fiuctuations, and even- 
tually annealed to zero. The route we followed (certainly 
not unique, or even the most efficient) goes through the 
mapping to a highly constrained Ising-like system, some- 
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what reminiscent of the Hopfield-Tank mapping of TSP 
as a neural network (See also Ref. 14.) Formally, 

each configuration of the system (a valid tour) is associ- 
ated to a N X N matrix T with 0/1 entries in the follow- 
ing way: For every directed tour (an ordered sequence of 
cities) Tij = 1 if the tour visits city i immediately after 
city j, and T^.j = otherwise. T has N entries equal 
to 1, all other elements being 0, and obeys the typical 
constraints of a permutation matrix, i.e.: a) All diagonal 
elements vanish, Ti^i = 0; b) There is a single 1 in each 
row i and in each column j. For the symmetric TSP 
problem we want to consider (a TSP where the distance 
matrix is symmetric dij = dji), a directed tour repre- 
sented by a T, and the reversed tour, represented by the 
transposed matrix T*, have exactly the same length. It is 
convenient, as will be apparent in a moment, to introduce 
the symmetric matrix U = T + as representative of 
undirected tours, the non-zero elements of U being given 
by Uij = 1 if i is connected to j in the tour, i being 
visited before or after j. One can be readily convinced 
that there is no loss of information in working with U 
instead of T. Given the matrix U there is no ambiguity 
in extracting the directed tour, represented by T, which 
originated it. The length of a tour can be expressed in 
terms of the Uij, as follows: 

Hpot {U) = ^^dij ili^j ^^d^jU^^j , ( 1 ) 

(y) (ij) 

where dij = dji is the distance between city i and city j, 
and {ij) signifies counting each link only once. Hpot is 
the potential energy we will finally seek to minimize. 

Since there is no natural physical kinetic energy in the 
problem we have to devise a suitable one. This choice is 
arbitrary, and many simple forms, such as transposition 
of two neighboring cities in a tour, could provide one 
or another kind of quantum fluctuations. Reasonably, 
however, the choice of Hkin should also encompass the 
important elementary "moves" of the problem, determin- 
ing which configurations are to become direct neighbors 



of a given configuration, a choice which in turn influences 
the problem's effective landscape [fsf . A very important 
move that is often used in heuristic TSP algorithms is 
the so-called 2- opt move, consisting in eliminating two 
links in the current tour, (ci — > C2) and (ci' — > C2')j 
and rebuilding a new tour in which the connections are 
given by (ci — > ci') and (c2 C2')- This is illustrated 
with a 8-city example in Fig. 1, were (left part) the tour 
(1^2^5^8^7^6^3^4), represented by the 
matrices Tin and [/jn, is rebuilt by eliminating the two 
links {ci — 2 ^ C2 = 5) and (ci' = 3 — > C2' = 4), 
and forming (see right part of Fig. 1) two new links 
(ci = 2 -> cv = 3) and (c2 = 5 ^ C2' = 4). The 
matrices Tfin and Ufin in the lower right part of Fig. 1 rep- 
resent the final (directed and undirected) tour obtained 
after the 2-opt move. A 2-opt move implies that a whole 
section of the original tour (between the two removed 
links (ci — > C2) and (cy — > C2') ) gets reversed in the new 
tour, yielding a long series of bit flips in the matrix T 
(dotted circles in Fig. 1). If we associate a spin variable 
+1 (—1) to each entry 1 (0), and represent a 2-opt move 
in terms of spin-flip operators acting on the configuration 
Tin, we would need a whole string of spin- flip operators 
(a non-local object) to enforce a trivial operation - re- 
versing a piece of tour - which does not affect the tour 
length at all. The advantage of working with U is that 
it represents, in a symmetric way, the direct and the re- 
verse tour, so that all the entries corresponding to the 
section of reversed tour are completely untouched. The 
whole 2-opt move, when working with U matrices, can 
be represented by just four spin-flip operators: 

where, by definition, each S^^ flips an Ising spin variable 
(defined as 'S'^-^-^ = {'iUij — 1) = ±1) at position (z,j) 
and at the symmetric position (j, i), i.e., = SfjS^^. 

The quantum Hamiltonian for the TSP which imple- 
ments the 2-opt moves is: 



Htsp = HpotiU) + Hi 



kin 



H.c. 



(2) 



(ij) {i'j'} 



Here the quantum coupling F is a real positive function 
depending, in principle, on the links, as well as, in the 
annealing problem, on time. The link dependence of F 
can be restricted to a dependence on the distances be- 
tween the cities involved in the two links created by the 
2-opt move, in such a way as to realize a neighborhood 



pruning |l2l |. by discouraging (or forbidding altogether) 
the creation of links between distant cities. In so doing, 
we mapped the symmetric TSP problem onto a highly 
constrained Ising spin problem with N{N — l)/2 sites - 
one for each pair of (j,j), with i > j, due to the sym- 
metry (i,j) (jj*), which we denote by -, with 
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FIG. 1: Left: Representation of an 8-city tour, with the 
corresponding matrix Tin and Uh-i = Tin + fii. Right: The 
final tour obtained when a 2-opt move is performed, with a 
whole section reversed (dotted line). The matrices Tfln and 
[/fin are shown, the circles indicating the entries that have 
been switched (0 ^ 1) by the 2-opt move. The dotted circles 
in Tfln are entries related to the trivial reversal of a section of 
the tour. 



a four-spin-flip kinetic term providing the 2-opt quan- 
tum fluctuations. The potential energy Hpot appear- 
ing in Htsp (tour length), represents, in the Ising lan- 
guage, a (random) external magnetic field at each site 
depending on the inter-city distance dij. The con- 
straints on the Hilbert space are such that the matrix 
IJi^j — {S^^ + l)/2 represents a valid tour. In particu- 
lar, the system lives in a subspace with a fixed magne- 
tization - exactly N spins are up {Uij = 1) among the 
N{N — l)/2 - and the 2-opt kinetic term conserves the 
magnetization. 

As in the past 0, we will not attempt an actual 
Schrodinger annealing evolution of the quantum Hamil- 
tonian proposed - out of the question due to the large 
Hilbert space. On the contrary, we shall address the 
quantum problem by Quantum Monte Carlo (QMC), 
where annealing will take place in the fictitious "time" 
represented by the number of MC steps. Path-integral 
Monte Carlo (PIMC) provides a natural tool, due to 
its simplicity. However, attacking Htsp by PIMC 
meets with a first difficulty, namely Trotter discretiza- 
tion of the path- integral [l6j. That requires calcula- 
tion of the matrix elements of the exponential operator 
(C'l exp — ei?fci„|C) between arbitrary configurations |C) 
and |C") of the system a complicated affair for the 
Hkin in Eq. To circumvent this difficulty without 

giving up the simplicity of PIMC, we introduce a dras- 
tic simplification to our kinetic energy term, replacing it 
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FIG. 2: Average residual excess length found after Monte 
Carlo annealing for a total time r (in MC steps), for the 
N = 1002 instance prl002 of the TSPLIB. Notice how quan- 
tum annealing (QA) provides residual excess lengths decaying 
faster than classical annealing (SA). 



altogether with a standard transverse Ising form: 



H' 



TSP 



s- 



{hi) 



H.c] , 



(3) 

This form is trivially Trotter-discretized, as in the stan- 
dard Ising system in transverse field since the spin- 
flip term acts independently on single spins at each time- 
slice. Strictly speaking this simplified form of kinetic 
energy would no longer fulfill the constraint to take a 
valid tour to another valid tour. However, so long as we 
constrain the MC dynamics strictly within the valid tour 
subspace - a restriction that comes automatically if we 
use exclusively 2-opt moves in the MC algorithm, and no 
single spin-fiip moves - that problem does not arise. That 
this way of bypassing the difficulty in treating the original 
Hamiltonian |(2Jl should eventually produce a working QA 
scheme is a priori far from obvious. We will eventually 
find that it indeed does. The simplified single spin-fiip 
kinetic term Hkin in Eq- ® enters only in calculating 
the weights of configurations of the Trotter discretized 
system, and not in the actual MC dynamics, which relies 
on the 2-opt moves and hence conserves the constraints. 
The details of the remaining implementation are identi- 
cal to those reported for the random Ising case 0, 0, 0| ■ 

For a direct test of our QA algorithm against SA we 
chose a standard benchmark TSP problem, namely the 
printed circuit board instance prl002 of the TSPLIB 
[3. It is a structured TSP problem with TV = 1002 
cities whose optimal tour length Lopt is known exactly 
by other ad-hoc algorithms 18] to be Lopt = 259045. 
Our implementation of SA was a standard Metropolis 
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MC with a temperature schedule starting from Tq and 
going hnearly to zero in a MC time r. We chose an op- 
timal initial temperature Tq by first performing several 
SA runs with different (short) cooling times r starting 
from sufficiently high temperatures (the ergodic region). 
This identifies an approximate "dynamical temperature" 
Tdyn below which the cooling curves for different t's start 
to differ. For prf002, we obtained Tdyn ~ 100. As it 
turns out, the optimal Tq for SA approximately coin- 
cides with Tdyn, Tq ~ Tdyn, an irnplementation feature 
known for TSP SA as cold staHs [ij. For QA, we im- 
plemented for the Hamiltonian Htsp a similar PIMC to 
that used previously|^ (ij at a fixed low temperature T 
(we used T = 10/3, see below). The quantum model is 
mapped onto a classical model with an extra imaginary- 
time dimension, consisting of P ferromagnetically cou- 
pled re plic as of the original spin problem, at temperature 
PT 0, [III (we used P = 30). Since QA requires initial 
configurations equilibrated at temperature PT , an ob- 
vious choice is to take PT ~ Tdyn, i-e., PT = 100 for the 
prl002 Finally, the transverse field F is annealed 

linearly in a MC time r from an initial value Fq = 300 to 
a final value of zero. In both SA and QA, we used exclu- 
sively 2-opt moves, with a static neighborhood pruning 
[l^ . which restricts the attempted 2-opt moves by al- 
lowing only a fixed number M (we used M = 20) of 
nearest neighbors of city j to be candidates for j'. Our 
MC step consisted of MN attempted 2-opt moves (for 
QA, in each of the P replicas). In both SA and QA, 
we averaged the best tour length found over up to 96 
independent searches. 

Fig. 2 shows the results obtained for the residual error 
(average best-tour excess length) upon annealing for a 
total MC time r, eexc = {Lbestir) - Lopt)/Lopt, both 
with SA (filled squares) and with QA (open circles). As a 
reference, the best out of 1000 runs of the Lin-Kernighan 
algorithrn (one of the standard local-search algorithms 
for TS P I12I I gives a percentage excess length of ~ 
1.57% [ia|7dashed fine in Fig. 2). The resuhs show that 
QA anneals more efficiently, reducing the error at a much 
steeper rate than SA. Moreover, even accounting for the 
extra factor P in the total CPU time required by QA 
(rightmost open circles), QA is still more convenient than 
SA at large r and small excess lengths. The similarity 
with the previous results on the random 2D Ising magnet 
is striking. We argued in Ref. that QA is faster 
for the random Ising case due to the superior ability of 
quantum physics to cross barriers through Landau-Zener 
tunneling, as compared to classical physics requiring for 
them to be overcome thermally: Such feature, although 
by no means a general property, is apparently shared by 
the TSP. The upward curvature of the data in Fig. 2 is 
likely to signal a logarithmically slow annealing for both 
SA and QA 0. Also worth mentioning is the effect of 
the finite value of P, which is likely to be responsible for 
a saturation effect of the QA data, as shown in Fig. 1 of 



Ref. la for the random Ising case. 

Current ad-hoc TSP algorithms do better than either 
annealings - SA and QA being generic tools - for the 
same CPU time. Moreover, it is known that combining 
SA with local search heuristics provides superior results 
to pure SA for the TSP problem Nevertheless, the 
absolute quality of QA and its success in a fair compari- 
son to SA strongly encourages to further applications of 
QA as a general purpose optimization technique, and to 
possible improvements of the bare Q A scheme by combin- 
ing it with other local heuristics, in the spirit of Ref. [20j . 
Equally instructive will be to experiment with other QA 
schemes, for instance Green's Function QMC, which are 
able to cope with the 2-opt Hkin constructed in Eq. (O or 
with other sources of quantum fluctuations. That broad- 
ening of the project we must leave for future studies. We 
believe that gaining further experience with the effects 
of artificially introduced quantum fluctuations in classi- 
cal complex problems represents a very promising and 
challenging route, particularly in view of future develop- 
ments in quantum computation. 
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